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VARIABLE  SUCCESSIVE  OVER -RELAXATION 

by 

Leland  Kitchin  McDowell 

Department  of  Mathematics 

University  of  Illinois,  1961 

This  thesis  investigates  the  solution  of  linear  systems  by 
extrapolated  Gauss-Seidel  iteration  using  a  multiplicity  of  extrapolation 
parameters,   Phis  is  a  generalization  of  the  method  of  successive  over- 
relaxation.,  which  uses  a  single  scalar  extrapolation  parameter.  The  linear 
systems  considered  are  those  which  arise  in  the  numerical  solution  of 
boundary  value  problems  for  self-adjoint,,  elliptic  partial  differential 
equations , 

In  Chapter  2  it  is  shown  that  the  use  of  two  appropriately 
chosen  scalar  extrapolation  factors  yields  an  iteration  having  a  higher 
rate  of  convergence  than  SOR,  and  formulas  are  derived  for  choosing  the 
factors  optimally,,  Also,  it  is  shown  that  by  the  use  of  extrapolation 
matrices,  an  iteration  can  be  constructed  whose  matrix  is  nilpotent,  i„e„5 
whose  rate  of  convergence  is  infinite. 

Chapter  3  considers  a  more  limited  class  of  linear  systems  for 
which  a  certain  sort  of  spectral  decomposition  is  possible.  For  the 
solution  of  such  systems  the  SE"  and  VSEJ-  methods  are  introduced  and  shown 
to  be  equivalent  to  several  simultaneous  extrapolated  lauss-Seidel  iterations 
on  certain  sabspaces^  each  with  a  different  extrapolation  factor  or  set  of 
factors o   Theoretical  and  experimental  results  are  presented  which  show  that 
SEIj  which  requires  less  work  per  iteration  than  SOR?  has  the  same 
asymptotic  rate  of  convergence  and.  for  most  starting  vectors,  an  improved 
actual  rate  of  convergence.   Experimental  evidence  is  presented  showing 
that  certain  versions  of  VSE7  have  a  higher  asymptotic  rate  of  convergence 
for  the  problems  considered  than  30R„ 


1.   INTRODUCTION 

1 „ 1  Relaxation  Methods 

The  use  of  relaxation  for  the  solution  of  linear  systems  is 
documented  as  early  as  1823;,  when  Gauss  [3]  wrote  favorably  of  his 
experience  with  such  a  procedure,,  which  he  called  indirect  elimination 0 
The  la cob i  Iteration  for  the  solution  of  a  linear  system 

AZ  --    K  (l.l.l) 

where  A  is  symmetric  was  discussed  in  1846  by  Jacobi  [7],  who  used  plane 
rotations  to  increase  the  diagonal  dominance  of  A  and  hence  to  improve 
convergence, 

In  1862,  Seidel  [10]  considered  the  method  now  known  as 
Gauss-Seidel  iteration  for  the  case  in  which  A  is  m  x  n  with  m  >  n„ 
Seidel  proved  that  his  iteration  applied  to 

AT  A  Z  =  AT  K  L.1.2) 

converges  to  a  point  which  best  satisfies  (1,1,1)  in  the  sense  of  least 
squares  and  hence  to  the  unique  solution  if  A  is  square  and  non-singular, 

Seidel  noted  that  if  there  exists  a  subsystem  of  k  unknowns 
and  k  equations  in  which  the  unknowns  within  the  subsystem  are  coupled  to 
only  a  few  outside  unknowns,  then  it  is  profitable  to  relax  repeatedly 
the  residuals  of  the  subsystem  until  they  become  small  while  treating  the 
outside  unknowns  as  constants,  Thus  he  anticipated  modern  block  relaxation 
methods^  in  which  a  number  of  residuals  are  relaxed  simultaneously, 


*  Numbers  in  square  brackets  refer  to  items  in  the  List  of  Referent 
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Both  the  Jacobi  and  the  Gauss-Seidel  iterations  are  known  to 
be  convergent  if,  for  example,,  the  coefficient  matrix  is  either  strictly 
or  jrreducibly  diagonally  dominant  [  5  ]  .- 

Experience  with  hand  computation  showed  that  it  is  often 
profitable  to  over-relax,  and  intuition  indicates  that  the  greatest 
advantage  lies  in  relaxing  the  largest  residual  at  each  stage.   If  this 
is  done,  then  the  judgment  of  the  relaxer  intervenes  to  alter  the 
algorithm  from  one  iteration  to  the  next,  making  the  iteration  non-cyclic. 

The  use  of  electronic  digital  computers  permits  the  application 
of  relaxation  methods  to  linear  systems  of  very  large  order,  but  such 
machines  are  best  suited  to  cyclic  iterative  methods,  in  which  the 
algorithm  once  defined  is  repeated  without  alteration  at  each  iteration,, 
For  linear  systems  arising  in  the  discretization  of  boundary  value 
problems  for  self-adjoint,  elliptic  partial  differential  equations, 
Frankel  [4]  and  Young  [12]  in  1950  described  such  a  cyclic  procedure 
for  extrapolated  Gauss-Seidel  iteration,  often  called  successive  over- 
relaxation  (SOR) o   Young  showed  that  for  certain  orderings  of  the 
unknowns,  which  he  termed  consistent  orderings,  a  functional  relationship 
exists  between  the  extrapolation  factor  w  and  the  rate  of  convergence  of 
SOR,   Using  this  relationship.  Young  derived  formulas  for  the  optimum 
value  of  w  and  for  the  rate  of  convergence  of  the  resulting  iteration, 
which  he  showed  to  be  substantially  faster  than  Gauss-Seidel  iteration, 
especially  when  the  latter  is  slow0 

Ybung; s  results,  which  applied  to  point  iterations,  were 
extended  to  block  iterations  in  195&  by  Arms,  Gates,  and  Zondek  [l]. 
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Few  results  are  known  concerning  the  use  of  extrapolation  in 
case  no  consistent  ordering  of  the  unknowns  exists,  but  Ostrowski  [9] 
proved  in  195^  that  if  A  is  Hermitian„  then  SOR  converges  if  and  only  if 
A  is  positive  definite  and  0  <  w  <  2„ 

In  1959  Golub  [6]  introduced  the  Chebyshev  semi-iterative 
method  and  the  modified  (or  cyclic)  Chebyshev  semi-iterative  method. 
The  latter  is  a  variation  of  SOR  asing  a  particular  consistent  ordering 
of  the  unknowns  and  a  sequence  of  extrapolation  factors,,  a  different 
pair  of  factors  being  used  for  each  iteration.   The  cyclic  Chebyshev 
semi-iterative  method  provides  an  improved  average  rate  of  convergence 
over  SOR,  although  the  asymptotic  rates  of  convergence  of  the  two 
iterations  are  the  same, 

1,2  Variable  Extrapolation  Factors 

Whereas  SOR  uses  the  same  scalar  extrapolation  factor  for  each 
unknown  of  the  linear  system  being  solved,  this  thesis  investigates 
extrapolated  Gauss -Seidel  iteration  using  extrapolation  factors  which 
vary  from  one  unknown  to  another  or  from  one  block  of  unknowns  to 
another.   The  use  of  matrices  as  extrapolation  parameters  is  also 
investigated.   We  call  such  iterations  variable  successive  over -relaxation 
(VSOR) 

The  linear  systems  considered  are  those  for  which  the  use  of  a 
single  extrapolation  factor  is  known  to  be  profitable ;  e0g„_,  linear 
systems  arising  from  the  discretization  of  boundary  value  problems  for 
self-adjoint  elliptic  partial  differential  equations 


k 

Chapter  2  begins  by  introducing  2-factor  VSOR,  an  extrapolated 
Gauss-Seidel  iteration  using  two  scalar  extrapolation  factors.   Formulas 
for  the  optimum  pair  of  factors  are  obtained,  and  it  is  shown  that  the 
rate  of  convergence  of  2-factor  VSOR  is  greater  than  that  of  SOR, 

In  Section  2  of  Chapter  2  the  concept  of  extrapolation  matrices 
is  introduced,  and  it  is  shown  that  by  their  use  an  extrapolated 
Gauss-Seidel  iteration  can  be  constructed  whose  iteration  matrix  is 
nilpotent,   Hence,  in  the  absence  of  rounding  errors,  the  exact  solution 
is  obtained  after  a  finite  number  of  iterations.   Moreover,  it  is  shown 
tnat  if  the  non-null  blocks  of  the  associated  block  Jacobi  matrix  are  all 
square  of  order  r  and  have  a  common  basis  of  eigenvectors,  then  VSOR  with 
extrapolation  matrices  is  equivalent  to  r  simultaneous  extrapolated  Gauss- 
Seidel  iterations  each  using  a  sequence  of  scalar  extrapolation  factors , 
Thus,  a  decomposition  theorem  for  VSOR  is  obtained. 

In  Chapter  3  we  consider  the  linear  system  resulting  from  the 
discretization  of  the  Dirichlet  problem  for  Poisson's  equation  on  a 
rectangle,   The  resulting  rq  x  rq  linear  system  is  the  model  problem  for 
whose  solution  the  sequential  extrapolated  implicit  method  (SET. )    is 
introduced.   This  method  accelerates  Gauss-Seidel  iteration  by  means  of  a 
shift  of  origin,  which  requires  less  work  per  iteration  than  extrapolation 
by  a  scalar,   It  is  shown,  however,  that  SEI  is  actually  equivalent  to  the 
use  of  certain  extrapolation  matrices,  and  the  decomposition  theorem  of 
Chapter  2  is  then  applied  to  show  that  SEI  is  equivalent  to  performing 
SOR  simultaneously  in  each  of  r  subspaces,  a  different  scalar  extrapolation 
factor  being  used  in  each  subspace.   It  is  then  shown  that  the  asymptotic 
rates  of  convergence  of  SEI  and  SOR  are  identical  when  the  optimum 
acceleration  parameter  is  used  for  each,  but  that  the  actual   rate  of 


5 
convergence  of  SEI  is  always  at  least  as  high  as  that  of  SOR  and  for 
certain  starting  vectors  is  substantially  higher, 

In  Section  3»3  cyclic  Chebyshev  SEI  is  introduced.  The 
relationship  of  this  iteration  to  SEI  is  analogous  to  the  relationship  of  the 
cyclic  Chebyshev  semi-iterative  method  to  SOP. : 

Section  3.^  presents  numerical  results  comparing  the  actual 
rates  of  convergence  using  various  starting  vectors  of  SOR,,  SEI ,  the 
Chebyshev  semi-iterative  method.,  and  cyclic  Chebyshev  SEI  , 

Finally,  the  variable  sequential  extrapolated  implicit  method 
( VSEI )  for  the  solution  of  the  model  problem  is  introduced.   Several 
versions  of  this  iteration  are  investigated,  each  of  which  uses  the  criteria 
of  Chapter  2  for  the  construction  of  nilpotent  iteration  matrices 
together  with  the  decomposition  theorem  to  make  the  VSEI  iteration  matrix 
nilpotent  on  certain  q-dimensional  subspaces  of  rq-space0   Experimental 
evidence  is  presented  which  shows  that  VSEI  is  asymptotically  faster  than 
SOR  for  the  model  problem- 


VARIABLE  SUCCESSIVE  OVER -RELAXATION  FOR  LINEAR  SYSTEMS 
WHOSE  COEFFICIENT  MATRICES  ARE  CONSISTENTLY  ORDERED  S-MATRICES 


2,1  Basic  Concepts 

Throughout  this  thesis  we  will  be  concerned  with  the  iterative 
solution  of  the  linear  system 


AZ  ^   K 


(2.1.1) 


where  A  is  a  matrix  usually  of  large  order  and  sparse;  i.e.,  having  only 
few  non-zero  entries.   A  detailed  development  of  the  ideas  outlined  in 
this  section  may  be  found  in  [11,  Chap.  1,3]  or  [13,  Chap.  1] . 

Our  results  will  apply  in  particular  in  case  A  is  an  S-matrix; 
i.e.,  A  satisfies 

(i)  A  is  real  and  symmetric 
'ill   a .  .  <  0  for  i  ^  j 

id  - 

(iii)  A  is  irreducible;  i.e.,  there  exists  no  permutation 
matrix  P  such  that 


T 
PAP  = 


A 


11 


0 


A21   A22 


where  the  diagonal  submatrices  are  square 
(iv)  A  is  diagonally  dominant;  i.e., 


li 


>  £ 
d# 


ij 


(v)  The  inequality  in  (iv)  is  strict  for  at  least  one  i 
S-matrices  arise,  for  example,  in  the  discretization  of  self-adjoint, 
elliptic  partial  differential  equations. 
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A     stationary,    linear   iteration  for  the   solution  of   (2„l.l)    Ls 
an  aff ine  transformation 

?m4l.HZm+F  (2.1.2) 

with  the  property  that 

Z  -  HZ  +  F  ['2.1.3) 

The  error  in  the  m-th  iterate  is  E  "  --  Z  -  Z   ,  and  subtraction  of  (2.1,3) 
from  (2.1.2)  shows  that 

m+1   TTZ*'  m   rTnn-JJf  o 


E =  HE   =:  H 


1fT°  1.10 


Hence;  the  convergence  properties  of  the  iteration  (2,1.2)  depend 

only  upon  the  matrix  H„   In  particular,   lim  E   =  0  for  arbitrary  E 

'   m-»<» 

m 
if  and  only  if  lim  H   is  the  null  matrix.   This  is  the  case  if  and  only 
m->°° 

if  all  eigenvalues  of  H  have  modulus  less  than  one.   The  spectral  radius 

of  H,  denoted  by  p(H),  is  the  maximum  of  the  moduli  of  the  eigenvalues  of 

H.  and  the  (asymptotic)  rate  of  convergence  of  the  iteration  (2,1,2)  is 

R  =  -In  p. 

A  partitioning  of  A  in  which  the  i-th  diagonal  submatrix  A 

is  r,  X  r.  induces  a  partitioning  Z  '   =  [Z  .  \  Z  _  J  — |Z    of  the  unknown 
11  ^  to  1  1   2  '   '   q' 

vector  Z  such  that  the  block  Z.  contains  r,  of  the  unknowns,   Tne  constant 

1  1 

vector  K  is  similarly  partitioned. 

If  A  is  an  S-matnx,  then  A  is  positive  definite,   Eact:  A,  .  is 

*  11 

therefore  also  positive  definite  and  hence  non-singular.   The  block 
Jacobi  iteration  corresponding  to  the  particular  partitioning  of  A  is 

Z  nifl  =  -  Z  A"1  A  .  Z  m+  A"1  K.  ;  1<  i  <  q,  (2.1.U) 

J+1   ii   lj    j     il   i     "   " 


where  the  partitioning  is  assumed  to  be  such  that  A. .  is  r.  x  r.. 

ij     1    J 

Throughout  this  thesis  M  will  denote  the  matrix  for  the  iteration  (2*1.^); 

i0e„ 3   M  denotes  the  block  Jacobi  matrix  derived  from  A, 

The  block  Gauss-Seidel  iteration  differs  from  the  block  Jacobi 

in  that  as  soon  as  a  block  Z  .    of  the  next  vector  iterate  is  calculated, 

1  ' 

it  and  not  Z  .  is  used  in  the  computation  of  later  blocks  Z  .  _  , 
l  l+l 

Z  . '   ,  etc.   Thus, 

14-2 

Zm+1  =   -    Z.   AT?"  A..  Zm+1   -    .X.  A"1  A  .  Z  m 

l       j<i  li   ij    j     j  >i   li   ij    j 


(2.1.5) 


+  A_1  K.  .    1  <  i  <  q. 


If  P  is  a  permutation  matrix,  then  the  coordinate 
transformation 

P  A  PT  Z  =  P  K  (2.1.6) 

0  0  0 

corresponds  to  a  reordering  of  the  unknowns.   The  ordering  o  is 

T 
consistent  if  A  =   P  A  P  is  block  tri-diagonal  (i.e.,  A. .  is  null  if 
ooo  '   ij 

j|>l)  with  square  diagonal  blocks,  and  if  so,  then  A  is  said  to  be 
consistently  ordered.   The  block  Jacobi  matrix  M  derived  from  A   is  said 
to  be  consistently  ordered  if  A  is.   Throughout  this  thesis  A  and  M 
will  be  assumed  consistently  ordered  already  and  the  subscript  o  will  be 
omitted, 

If  M  is  consistently  ordered,  it  is  of  the  form 


q-lq-1  q-lq 


a"1  A   _    0 
qq  q^q-i 


o    B 


q-lq 


so  B.  .  =  0  for  ji-j|  =f=  1°   Kence  M  =  L  +  U;  where  L  is  strictly  lower 
triangular  and  U  is  strictly  upper  triangular „   Moveover,,  if  M  is 
consistently  ordered  and  |i  is  an  eigenvalue,  then  -u  is  an  eigenvalue 
of  the  same  multiplicity. 
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If 


D  = 


-1  -1/2    -1/2 

then  M  =  I  -  D   A.   Hence  M  is  similar  to  I  -  D  '   AD  '  ,    and  it 

follows  that  this  is  symmetric  because  the  square  root  of  a  positive 

definite  symmetric  matrix  is  symmetric .   Hence  all  of  the  eigenvalues 

of  M  are  real, 

2 , 2  Successive  Over-Relaxation  (SOR) 

Define  the  intermediate  value 


_  m+-         _   . 

Z  .  2  =  B.  .  _  Z  m+]   +  B.  . 

1      i,i-l   l-l    1,1+1   l+l    11   1 


z  m  .  +  a:1  k 


(2.2.1) 


and  extrapolate  by  a  fixed  scalar  W; 


1 

Z  m+1  =  "[Z  ^  -  Z  m]  +  Z  ?,   1  <  i  <  q 

i         j_       i      i 


(2.2.2) 


Then 


Z   m+1  -  u>[B.     -    ,    Z  ^   +  B.     .    .    Z    . 
1  i,i-l        l-l  i,i+l        l+l 


H  WAT.  K.,   1  <  i  <  q< 
11   1     —   - 


m  ]  +  (i-«)  z  m 


(2.2.3) 


The  iteration  (2.2.3)  is  the  successive  over-relaxation  method  with 
extrapolation  factor  <£.   The  matrix  for  the  iteration  (2.2.3)  is 


11 


£     =»   [I   -  coL]"1    [uU  +    (l    -   co)l] 

CO 


Let  0  <  Li  <  i-t  < <  Li  =  (J.  be  the  positive  eigenvalues  of 

M,  and  let  p  be  the  multiplicity  of  zero  as  an  eigenvalue.  For  1  <  i  <  ly 
let 


CP.(\)  =  \      -    [2(1-W)  +  CO   Li  ]  X   +    (l-w)   . 
Then  tne  characteristic  polynomial  of  £  is 

CO 

$;\)  ^  [(i  -  co)  -  \]*   n  cp.(\)  . 

isl   X 

The  maximum  of  the  moduli  of  the  roots  of  cp  (A.)  is  minimized  as 

JO 

a  function  of  lo  when  the  roots  are  real  and  equal;  that  is,,  when  the 
discriminant  of  9/A)  is  zero„   This  is  the  case  when 


co  =  co  s 
b 


+  \|l  -  |I 


For  co  -  co   the  modulus  of  the  roots  is 
b7 


pSOR 


1  -   1-L  -  Li 


i  +  \fT7^ 


If  the  discriminant  of  cp  (\)    is  negative  or  zero  for  a 
particular  value  of  co,,  then  the  discriminant  of  cp.(A,).,  1  <  i  <  l0    is 
also  negative  or  zero.   It  follows  that  for  co  ^  co  the  roots  of  <E>(A.)  all 
lie  on  a  circle  of  radius  pariDp  and  hence 

min  pUJ  ^  p(j^  )  =pS0R  . 
co  b 

This  completes  the  summary  of  known  results. 
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2.3  VSOR  With  Two  Extrapolation  Factors 

In  this  section  we  consider  a  linear,  stationary  iteration  using 
two  distinct  extrapolation  factors .   It  will  be  shown  that  for  certain 
consistently  ordered  S-matrices  an  optimally  chosen  pair  of  factors  used 
in  alternation  on  consecutive  blocks  yield  an  iteration  having  a  higher 
rate  of  convergence  than  SOR. 

For  any  two  non-zero  scalars  w  and  w  consider  the  iteration 


T1"1  =  u>  [B   ,Z^  +  B      »  m 
1      l   1,1-1   i-l    l 

(2.3.1) 


Z.1+X  =  u.[B.  .  ,  Z  n   +  B.  .  n  Z  m     I  +  (1  -  U>  )  Z  m 

-.,1+1   r+lJ    s      i'    l 


+  U).  A?1  K.,  1  <  i  <  q, 
l   li   i'    -   -  ^ 


where 


w        is  odd 

CO. 

1 w   if  i  is  even 


Let  Ju      denote  the  iteration  matrix  for  (2.3. l),  and  suppose  that  X 

1'  2       _ 
partitioned  like  Z  is  an  eigenvector  of  £   w  belonging  to  the  eigenvalue 

1'  2 
It  follows  from  (2.3. l)  that 

AX.  =  w.[\B.  .  n  X.  n  +  B.  .  _  X.  J  +  (1  -  w. )  X. ,  1  <  i  <  q, 
l     l    i,i-l   i-l    i,i+l  l+l  l   l'    -   —  ' 

or,  after  dividing  by  w. , 

(l  -  w. )  -  A 

0  -  [X.B.  .    X    +  B.      X.   1  +  [ —1 ]  X.  . 

1,1-1   i-l     1,1+1   l+l  w.  i 

'  '  l 

— »  \^  — "   — ■> 

This  linear  transformation  on  X  we  write  as  G  X  =  0.   If 


7,  = 


(1  -  u.)  -  \ 


i       w. 

l 


then 
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G 


q-i;q 


\B    .    y      I 

q;q-l   7q  q 


Here  I   is  the  identity  matrix  of  the  same  order  as  A..„ 
1  11 

\s    — *     — * 

The  equation  G  X  =  0  has  a  non-trivial  solution  if  and  only  if 

det  (Gj  =  0,  and  since  det  (G.)  is  a  polynomial  of  degree  n  in  A.;  its 

roots  are  the  eigenvalues  of  £     „  Assuming  A.  f   0  (an  assumption  which 

V    2 
will  later  be  justified;,  we  premultipy  G  by 


V 


1 
.  "2 


and  postmultiply  by 


A  -- 
r 


\    * 


X 


lit 

This  transformation  multiplies  each  subdiagonal  block  by  A.   ,  each 

1 

"2 
diagonal  block  by  \     }    and  each  superdiagonal  block  by  1.   Hence, 

1 

~2  n/ 
letting  7.  =  A.   7  . ,  we  have 


G  =  A.   G  A 


q-i,q 


B        7   I 
q,q-l     q   q 


But  det  (A.)  det  (A  )  =  X.   ,  and  so  the  following  lemma  has  been  proved, 


Lemma  2,3^1:   X.  ./•  0  is  an  eigenvalue  of  Jo      if  and  only  if  it  is  a 

1-  2 
zero  of  det  (G). 

Since  A. .  is  of  order  r.  the  vector  Z.  contains  r 
11  111 

components,  and  it  follows  that  w  is  used  with  r  =   Z   r.  components 

i  odd 

and  w  is  used  with  r  =    Z   r.  components. 
2  e    .       1 

1  even 

Lemma  2.3 >2;   Let  the  matrix  G'  be  obtained  by  replacing  each  diagonal 

\  1/2 
element  7.  of  G  by  (7  7  J    .   Then 


r   -  r 
o    e 


r   -  r 
e    o 


det(G)  =  7 


1 


det(G'  ) 
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Proof 


Let 


1 

yl      1 


r  = 

r 


1 
r2     2 


72I 


and 


r 


i 

r"2I 
72   1 


1 

7_2I 
rl       2 


1 
7_2I 


\ 


Then  G  =  T0   G"  r  -  and  so  det  (g)  =  det  (r.)  det  (G' )  det  (r  ).     But 


det(r^)   det(rr)   -  71 


r      -   r  r      -   r 

o  e  e  o 

2  2 

7  ,    -which  proves   the   lemma 
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Lemma  2,3^3;   Let  U  ,    u.   ,    ..,,    u  be  the  positive  eigenvalues  of  M,  and 
let  p  denote  the  multiplicity  of  zero  as  an  eigenvalue.   For  1  <  i  <  £., 
define 

<P   (X)  =  X  -  [(1  -  wi)  +  (1  -  o>2)  +  ^uy^]  \  +  (i  _  Wi)(i  _  w2) 

Then  the  zeroes  of  det  (G* )  are  precisely  the  zeroes  of 


£ 


E  E 

[(1  -  ^)  -  ^]2  [(1  -  "2)  -  ^f 

i=l   i 


Proof : 


G'  = 


l/2   ]/2 
is,  except  for  the  diagonal  blocks  y   '      y   '      l.}    identical  to  M.   Hence 

1/2   1/2 

det(G!  )  =  0  if  and  only  if  7    7    =  \i.,    where  \x.    is  an  eigenvalue  of  M, 

l/?  l/? 
Moreover,  the  multiplicity  of  (7    7    -  u.)  as  a  factor  of  det  (G' )  is 

equal  to  the  multiplicity  of  u.  as  an  eigenvalue  of  M,   If  |J..  /  0,  then 

/  1/2  1/2     N  ,      N 

(7,    7„   +  M.J  is  also  a  factor  of  det  (G' )    and  therefore  so  is 
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2                2  2 

17!  yk  But  ^  72  -       --— 

from  wh  •■■,s  that  corresponding  to  the  21   non-  zero  eigenvali 

t 
of  M  are  2#  zeroes  of  di  I       ,^mei>  the  21   roots  of  II 

The  remaining  zeroes  of  det     are  the  ze: 

1  1        f, 


genvalues 


•  2  2     p 
r,  /0  -  01 


l-w]  I  -  m[  l-a>2)  -  M 


uj   a)   A. 


from  whicn  the  conclusion  of  the  lemma  follows „ 

p  +  (r  -  r  p  +  (r  -  r  ) 

o    e  e    o 

Lemma  2  -  3 and  p  a —  are  both  non- 
negative  integers. 


Proof'   r  -  r  =  n  =  p  +  2£„     Hence,   r  +  r  )  has  the  same  parity  as 
o    e  '  o    e ' 

p,  and  so,  therefore,  does  r   -  r    Then  p  +  (r   -  r   and 

o    e  e 

-  r  i  are  both  even,  and  so  p^  and  p^  are  integers, 
■  o    e'  12 

By  lemmas  2  3*2  and  2 „ 3 -  3 -  the  zeroes  of  det   5)  are  precisely 

Px  P2   I 

the  zeroes  of  [{l  -  w  )  -  Vj   [(i  -  w     K]  n  cp  '  \  •  ,   rh<    are 

exactly  n  such  zeroes  because  each  one  is  an  eigenvalue  of  £      by 

2 
lemma  2       Hence  p  <  p  and  p  <  p„   But  p  +  p    p.  so  both  p  and  p 

are  non-negaiiv- 


orera_2      The  cnaracteristic  polynomial  of  £,      is 

2 

P]  ^ 

$..■   ~  [(l  -       II  cp 

i 


f:      B;y    lemmas   2  2   3  2-    and   2  of  $>(> )    are 

:isely    the   eigenvalues  of  £.      ,  But  by  lemma     2-3-'4     <£(.  V1    is  a 

2 

nomialj  which  proves  the  theorem. 
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In  practical  problems  the  eigenvalues  of  M  are  usually  not 
Known,  but  an  interval  can  be  determined  in  which  the  eigenvalues  are 
known  to  lie.   We  therefore  consider  choosing  to  and  to  to  minimize  the 
maximum  of  the  moduli  of  the  roots  of 

cp  (\)  =  \2  -    [(1  -  u^)  +  (l  -  u>2)  +  ^2[i2]\   +  (1  -  wi)(l  -  w  ) 

where   it   is  known  only  that   0<u_<u-<u<l8 

The   discriminant  of  cp   (X)    is 

u 


D(n)  -  w1w2ki2[wiw2H2  +  2(1  -  w  )  +  2(1 


to  )  1  +  (to 
2;j  i-  v  2 


V 


Lemma 


2,3.5:    If  ^  >  l,  to  >  i,  D(u_)  <  0,  and  D(u)  <  0,  then  D(u)  < 


for  \i  <   |i  <  \i, 


Proof:   This  is  evident  from  the  graph  of  D(u). 


D(n) 


J*»  H 


Figure  1 


note  the  maximum  of  the  moduli  of  the  roots  of 

Then    max        ITS  minimal  if  and  only  if  to  and  to  are  such 
H  J  1      2 

that  D(u/  =  D(ff    0: 


Proof   The  equation  cp  (\)  =  0  can  be  written 


-  U  -  (1  -  u>  )][\  -  (1  -  u)  )]  =  u2  \ 
12  X 


Completing  the  square  oh  the  left  hand  side  yields 

(i  -  k)  -  (i  -  <*>  tp 


CO    to 

1?L 


\ 


(to      -   to   V 
2    J  2  l' 
a   U       A.    + 


i|CO    U) 

12 


Some  algebraic  manipulation  converts  this  to 


\pl 


i_  (X  .  1)  +  i 


u 


\  co      \  CO 
\     2     N  .1 


LO 


2     1 
=  u   \.  +  £ 


=1  2 


CO      \  CO 

o       \    i 


It  follows  that  the  roots  of  9  'M  =  0  are  the  abscissas  of  the  points  of 
intersection  in  the  (a.  a)  plane  of  the  parabola 


2    2,      2.1 


to 


, 


CO 


\  Cu  cu 

\|  1  2 


to 


■ 
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Figure  2 


(jj 


2  2 

77-  =  I,  the  vertices  of  the  parabolas  a ,  (p.)  are  at  the  origin.   For 


Ui 


1 


U) 


0) 


7^  1,  the  parabolas  intersect  on  the  a-axis,  the  vertices  lie  on  the 


1 


negative  \-axis,  and  the  abscissas  of  the  vertices  increase  toward  zero 

2 
as  u  increases,   For  fixed  —  >  1  and  ^-.^p  sufficiently  small,   the  line 

2  ~~ 
a     intersects  the  parabola   a  (u)  in  two  distinct  points.   As  w-,wp 

increases,  the  slope  of  the  line  decreases,  and  so  p—  decreases  until  the 

line  becomes  tangent  to  the  parabola.  At  this  point  the  roots  of  cp— (\) 


are  real  and  equal  to  \|  ( 1  -  w  )(l  -  w  ).   If  w  >  1  and  w  >  1  (as  will 
be  shown  to  be  the  case  of  interest),  then  further  increase  of  w  u 
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causes  the  roots  of  cpH\)  to  become  complex  and  to  increase  in  modulus c 


u 

2 

It  follows  then.,  tnat  for  —  fixed,  the  value  of  ^-.^  which  minimizes  p— 

w,       '  12  u. 

o  .— 

is  the  value  such  that  the  line  g  is  tangent  to  the  parabola  a,  (n)<> 

co 
2 
For  —  >  1  but  sufficiently  small,  the  line  a„   does  not  intersect 

the  parabola  a   (u_)  when- it  is  tangent  to  the  parabola  a,(|^)j  and  so  the 

CO 

roots  of  <?..(*■ )  are  complex  of  modulus  Wl  -  w,)(l  -  w0)o  As  — —  increases 


2, 


the  parabolas  a-,  [n)  and  a,  (u)  move  to  the  left  and  the  distance  between 


their  vertices  increases, 

CO  CO       CO 

2  2     2\ 

For  -—  sufficiently  large,  say  —  =  (— v,  ,  there  exists  a  value 

CO  •'W       CO   h 

1  11 

of  co  w  -which  makes  the  line  n  simultaneously  tangent  to  the  parabolas 

2  -        2 / ■— 

7     and  a,  >n).  Hence,  D(n)  =  D(^)  --■   0,  and  the  roots  of  <p—(?0  are  both 
I  *"~       1 


equal  to  U,l  -  w-.)(l  -  w„)  while  the  roots  of  cp   \;  are  botn  equal  to 

-  \|  1  -  w  )(l  -  w  )  „   Moreover.,  it  follows  from  lemma  2,3-5  that  the  roots 

of  cp  (X)  for  \x  <  [i   <  |i  are  complex  of  modulus  \|'l  -  w-.)(l  -  ^p 

cu 
2 

..creases  further  and  w,^  is  adjusted  so  that  the  line 
"J  12 

p , 
remains  tangent  to  the  parabola  a-,'-.';  the  point  of  tangency  moves  to  the 

left  because  the  parabola  moves  to  the  left0  Hence  min  p..  regarded  as  a 

function  of  U)-,w0  increases,  and  the  line  ceases  to  intersect  the  parabola 

2-  - 
a,  ,- ■) ,      It  follows  that   max     p|p  is  minimal  when  to  and  w  are  such 

b.  _  ~  : 

that  the  line  is  simultaneously  tangent  to  both  parabolas ,   which  is  the 

case  if  and  only  if  D  0, 
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o 
This  argument  shows  that  if  —  >  1,  then  the  pair  (oj  w  )  is 

1 

optimal  if  and  only  if  D(u)  =  D(u_)  =  0.   This  latter  condition  also 

00 

2 
characterizes  an  optimal  pair  such  that  —  <  1,  as  can  be  shown  by 


OJ 


1 


letting  —  decrease  from  1  in  the  above  argument.   Indeed,  since  cp  (X.)  is 


tc 


a  symmetric  function  of  w  and  co   it  follows  that  (oo  oo  )  is  an  optimal 


pa 


ir   if  and  only  if   (w     oo   )    is   alsc 


2 
Theorem  2C3.3:    Let  "     - , 

b        1+   u_u  +  \J(l-u_2)(l-JI2) 


w 


1   -   u  u  +  \|(l-u.2)(l-Ti2) 


,      and       p 


W    -   \ 


1-J2 


2S0R 


\Jl^?     +    \|l  -  ? 


Then 


pu  =  P2S0R  ^ 


max       _  p 


u_<u<^  w]_>u2    |  li  <  ^  <  I1 


M- 


00        =  OJ  j  00  =     00 

1   and  only      '  either <  -     >  or  <     ,  , ,      )     holds. 

»  OJ   =  00   /  \  0J  =  00 

2  b  I       2  -b 


S0R 


, )   If  w  denotes  the  optimum  extrapolation  factor  for  SOR  and 
p(£  ),  thenl  <  to  <  w  <  ST  and  pos„  <  p 


U> 


b  -     b  -     b 


"2S0R   -  MS0R' 


(iii)      uob   =  oob  =  oob   and  p^  =  pg0R  if  and  only    if  u_  =   0, 


Proof: 


(i)   By  Theorem  2,3.2,    max   _  p  is  minimal  if  and  only  if 


u_  <  u  <  u 


Dl  ,,  )  -  D(u)  =  0.   Since 


D(|i)   -  CdjUyi      [oo^u2  +   2(1-^)    +  2(l-w2)]   +    (oj^oj^ 


2  2-2 

is   quadratic   in  u-     with  roots  u_     and  |i   }   we   have 


(cu        I     U)        -         ) 

1 


CO      to 
1      2 


(w  -  to   V 

£     u     =        2     p 

CO  CO 

1  2 


(2.3.3) 


CO        -     CO 

Adding  plus   and  minus   the   positive   square   root  Li  u_  -•  — — — —     of   (2,3  =  3) 


to    [2  3-2'   yields 


CO        CO 

1      2 


(u  +  u_)     co^   -   Uco2  +   h 


and 


-  u_)     co^2  -  to^  +  h  =  0, 


whose   solutions   are 


and 


CO       ^    CO 

1       -b 


CO       =    CO 

2  b 


1 


"\ 


1     t      J    J     -    \|(l     -    U     )(l    -    ji     ) 


U) 


fc»    f- 


\J(1      -      ^)(1      -      M2)j 


(l-<*>  )(l-u>  )|  is  smaller  for  the  solution 


2h 


CO   =  CO 

1   -b 


CO   =  CO 

2    b 


An  analogous  development  using  the  negative  square  root  (i  (i  = 


of  (2.3.3)  yields  the  soluti 


on 


co       =    CO 

1  b 


CO       =    CO 

2        ~b 


Hence, 


min   /  max     _     p       >  =  ^(l   -  w   )  (l   -   cq   )    } 


CO        -     CO 

_1 2 

CO    CO 

1   2 


ind  some  algebra  shows  that  ^(l  -  ^h)(l  -  w,  )  =  p 


2S0R 


(ii)  Since  0  <  kL  <  n  <  1,  ™e  can  set  u_  =  sin  a  and  \x   =  sin  p, 


re  0  ^  a  <  p  <  -    Then 


CO    3 


-b  "'  1  +  sin  a  sin  P  +  cos  a  cos  P  "   1  +  cos  (p  -  a) 


co     = 


b        1   -   sin  a  sin  P  +   cos  OL  cos  P  "    1  +   cos    (p  +  a) 


and 


co     = 
b 


A 


1   -  \xc  1  +   cos  P 


Sine         C  a  <  0  <   £,    it  follows   that  cos    (p  +  cc)   ;_  cop  p  <  cos    (p   -  a), 

■ 

and  so   1  <  w     <  co     <  co   . 

b  -     b  ~     b 

If  u_  =  0,    then 


J2S0R 


i  -Ni  -  H2 


■N[ 


■  ? 


pS0R 


But   C—    p  <  0  for  u.  >  0,    and  therefore    *  <  p  for  u  >  0. 

riu_   MiS0R  U  ^2S0R  -  pS0R  L  - 
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(iii)  This  follows  at  once  from  the  formulas  which  define  the 

quantities  involved. 

It  is  of  interest  to  trace  the  movement  in  the  complex  plane  of 

the  roots  of  <P  (\)  for  fixed  w  end  w  satisfying  1  <  w  <  u  as  n 

decreases  from  +  «  to  zero.  It  follows  from  (2.3.2)  and  (2.3.3)  in  the 

proof  of  the  preceding  theorem  that  each  such  pair  (w  u)  )  corresponds  to 

some  pair  (u_>  u)  satisfying  0  <  u^  <  u. 


(l-"2) 


\J(l-u)i)(l-w2) 


Figure  3 


For  u  =  +  c»5  the  roots  of  cp  (k)   are  zero  and  +  <»„  As  j_  decreases.,  the 
roots  approach  each  other  until  for  u  =  u,,  they  are  equal  to 
Yl  -  un)(l  -  ^p^0  ^E  ^  decreases  from  u  to  u_,  the  roots  move  around  the 
circle  of  radius  m«,  1  -  w,  )(l  -  °0  in  the  complex  plane,  becoming  equal 
to  -  y'l  -  w)(l  -  w  )  when  u  =  ji.  As  u  decreased  from  u_  to  zero,  the 
roots  move  apart  on  the  real  axis  so  that  for  u  ■--.   0;  one  root  is  (l  -  u  ) 
and  the  other  is  (l  -  w  )„ 
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Theorem  2,3°^  Using  the  notation  of  lemmas  2„3ol  and  2„3.>3>  let 


tL  = 


u_   if  p  =  r  -  r 
1         '  o    e 


0   if  p  >  r  -  r 

'  o    e 


> 


Then 


0^   ,u>    }   =P2S0R  =       mln        pUco  ,W    }    • 

Proof:      If  p  =    |  r      -  r   |,    then  by  theorem  2<,3°1  the   characteristic 

polynomial  of  £       w     is 
1'    2 

P     ^ 
$(\)  =   [(1  -  w  )   -  M       n     cpn   (M    o 

^  i=l     •*! 

By  theorem  2„3o3?  the  maximum  of  the  moduli  of  the  roots  of  II  cp  (\) 

i=l  **I 

is  minimal  for  either  of  the  two  choices  for  (u)  u  )  specified  in  the 

theorem,  and  in  fact  the  roots  lie  on  a  circle  of  radius  Poono°  But  the 

2oUn 

remaining  roots  of  0(\)  lie  inside  this  circle  because 


|1  "  ^ul  <  Y-l  "  ->)(-'-  "  ^v)  a  P0nnRj  an^  the  conclusion  of  the  theorem 

follows  „ 

If  p  >  |  r  -  r  I  ,  then  both  (l  -  w  )  and  (l  -  w  )  are  roots  of 
1  o    e1  '  -b  b 

$(X),  But  these  are  the  roots  of  cp  (.\)  if  [i   =  0,  and  so  it  follows  that 

o'£,   w  )  is  minimized  by  taking  ^  =  0„ 
1-"  2 
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2 , h     General  VSOR :   Extrapo lation  Matrices 

Let  A  be  partitioned  so  that  each  diagonal  submatrix  A  .  is 

11 

square  of  order  r. ,  All  of  the  iterations  considered  in  this  thesis  for 
the  solution  of  (2.1.l)  are  of  the  form 


r  m+1 


=  ft.[B.  .   Z 


~  m+l 


+  B 


m 


i   i,i-l  i-1    i.i+1  i+-l 


]  4   (I  -  fi.)Z 


rn 


+  fl  A   K  ,  1  <  i  <  q 
l  li  l 


(2.U.1) 


where  each  ft.  is  a  non-singular  square  matrix  of  order  r.0   For  example, 

in  the  case  of  SOR  or  of  the  2  factor  VSOR  of  Section  2.3,  each  ft.  is 

simply  a  scalar  multiple  of  the  identity  matrix.   The  iterations  to  be 

investigated  in  Chapter  3  are  also  special  cases  of  (2„U,l),  although 

their  implementation  does  not  involve  explicit  multiplication  by  the  ft.'s 

i 

If 


o" 


o 


then  the  iteration  matrix  for  (2.U„l)  is 


r,  »  [i  -  nL]"  [nu  +  (i  -  n) 


where  L  +  U  -  M,  the  block.  Jacobi  matrix  derived  from  A,   If  each  ft.  is 

'  l 


non-singular  and   if  Z        -   Z.    the   true  solution  of   (2.1.l)5    then  Z 
also, 


rm+i  a  - 
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In  case  each  fi.  is  a  full  matrix,,  the  implementation  of  (2„40l) 
evidently  entails  considerably  more  work  than  extrapolation  by  a  scalar „ 
Suppose,  for  example,  that  r.  ~r,  l<i<q0   Z.    can  be  computed  most 
economically  (even  if  fi.  is  a  scalar  multiple  of  the  identity)  by  first 
computing 


1 

I™''  B.  .  ^4  B.    2*   +A;!if.  (2.U.2) 

l        i.i-l  l-l    i.i-ii-l    hi 


and  then  extrapolating: 

1 

2*1-  a.    [t*2  -   zmj  4  z  m  .  (2.U.3) 

11111  ' 

2  ""* 

;  2,4,3)  involves  r  multiplications  and  r(r+l)  additions,  and  since  Z. 

contains  r  unknowns,  it  follows  that  the  extrapolation  procedure  for 

2„^ol)  requires  r  multiplications  and  r+1  additions  per  unknown  of  the 

linear  system  (2.1„l)o   Extrapolation  by  a  scalar,  however,,  requires  only 

one  multiplication  and  two  additions  per  unknown.   The  direct  use, 

tnerefore„  of  full  extrapolation  matrices  is  computationally  advantageous 

only  if  a  large  increase  in  the  rate  of  convergence  can  be  obtained  over 

iterations  using  scalar  extrapolation  factors. 

Such  a  large  increase  is  indeed  possible.   In  fact,,  for  a 

properly  chosen  set  [fi.]  of  extrapolation  matrices,  £,  is  nilpotent  (i„e0! 

p'Ju  '  ='•  Oj  so  that  the  rate  of  convergence  of  (2.4,l)  is  infinite,   The 

following  theorem  gives  criteria  for  choosing  such  a  set  (fi„ }. 
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Theorem  2;U, 1,    jV  is  nilpotent  if  ft  ,    ft  , ft  satisfy  one  of 

(i)   fi.  -  I  and  ft.  =.  [I  -  B.  .   0.   B  .  J"1.  2  <  i  <  q 
1         i        i,i-l  i-I  i-l,i   J   -   —  ^ 

(ii  i     fi      =   I   and  ft.    =.    [I   -  B.    .    _fi.    _B.    .    .  ]      ,    1  <  i  <  q-1 
q  i  i,i+l  l+l   i+l,i        }        -       -  ^ 

(iii)      ft     =  I,    ft     =»  I,    and   for   some  m  such  that   1  <  m  <  q; 

ft.    =    [I    -   B.     .    .ft.    nB.    .     . ]"ls    2  <  i   <  m, 
l  i, i-I   i-I   i-I, i  -  ' 

ft.    =    [I    -   B.     .      fl,      B.    .     .J"1,    m+1  <  i  <  q- 
1  l,14l    1+1   i+l,i         '  -        -  ^ 


ft     -•    [I    -   B  ,ft      _B 


m 


B  ft        B 


-1 


m,m-l  m-1  m-l,m         m,m+l  m+1  nH-l,m 


Proof ;   Suppose  that 


x-»ilxI|--|x^ 


is  an  eigenvector  of  £,„  belonging  to  the  eigenvalue  K0      It  follows  from 

(2.4,1)  that  for  1<  i  <  q,  KK.    ..•  ft,  UB.    nX.  _  +  B.  .  nX.  ,  ]  4   I-ft  >X. 

-  y        l  i   1,1-1  i-I    i,i+l  i+1        i   i 

and  hence  JX  --■  0,  where  G  is 


-u  )-\        ft  B 

'll'    1J      1  12 

v 


2  21 


\ft   B 

q  q,q-i 


ft   ,B  . 

q-i  q-l,q 


q  q' 


Here  I.  is  the  identity  matrix  of  order  r  , 
i  i 
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There  exists  a  non-trivial  solution  of  GX  -  0  if  and  only  if 
det(G)  =  0,  and  since  detiG)  is  a  polynomial  of  degree  n  in  \,  it  follows 
that  its  roots  are  precisely  the  eigenvalues  of  £  „   We  now  show  that 
if  the  ft.'s  satisfy  (i),   ii;,  or  (iii),  then  det(G)  ~  +\    0 

If  (i)  is  satisfied,,  then  the  block  in  the  (l,  l)  position  of  G 

is  -VI  „   We  now  perform  the  block  elementary  row  operation  of  adding 

Q  13   times  the  first  rn  rows  of  G  to  the  next  r  rows ,   That  is ,  we  add 
2  21  1  2  ' 

n2B21(-M1)  to  Mi2B21  and  ^21^12   to  [(i  -  I^)  -  Xl^].      This  operation 
leaves  det' G;  unaltered  since  it  corresponds  to  premultipli cation  of  G 
by 


I 


21 


L 


o 


whose  determinant  is  obviously  1„ 

block  in  the  (2,l)  position  is  now  null,,  while  the  block 
in  the  (2,2)  position  is  [ft  B  ft  B   +  (i-fl  )]  -  \1  „   The  quantity  in 
square  brackets  vanishes  because  ft  satisfies  the  recursion  relation 
and  so  the  diagonal  block  in  the  .2,2)  position  is  -XI     Hence.,  we 
repeat  the  procedure,  adding  ^B   times  the  second  r  rows  to  the  next 
r  rows,   The  sub -diagonal  block  in  the  (3,2)  position  vanishes^  while 
the  diagonal  block  in  the  (3,3)  position  becomes  -M   since  ft  satisfies 

t  ■  .   P}  continuing  in  this  waj. ,  all  sub-diagonal  blocks  of  G  are  made  to 
vanish  (i.e.,  G  is  made  block  upper  triangular},  and  the  i-t'h  diagonal 
block  becomes  -\I.„   This  procedure  leaves  det(G)  unaltered,,  and 
det  (j)    t  \    „   Hence  £,  is  nilpctent. 
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case  the  ft.'s  satisfy  (ii),  we  proceed  by  block  column 

operations  starting  from  the  right  to  reduce  1  to  the  same  upper 

triangular  form. 

If  (iii)  is  satisfied,  we  proceed  by  block  row  operations 

starting  from  above  to  eliminate  the  sub-diagonal  blocks  in  positions 

(2,lV  (3,2),  (m,m-l)  and  by  block  column  operations  starting  from 

the  right  to  eliminate  the  remaining  sub-diagonal  blocks,   Again  the  i-th 

diagonal  block  becomes  -XI. .  Hence,  if  one  of  (1     Li)   or  ^iii)  is 

satisfied,  det(G)  -  +  A.  ,  and  so  JlL  is  nilpotent  in  all  cases, 

QoE0D„ 

If  Jo  is  nilpotent,  then  the  error  in  the  successive  iterates 

eventually  becomes  zero  in  the  absence  of  rounding  error,  and  so  the 

exact  solution  to  (2al„l)  is  obtained  after  a  finite  number  of 

iterations.   By  the  Cay ley -Hamilton  theorem  at  most  n  iterations  are 

required,  and  in  fact  it  can  be  shown  tnat  q  iterations  suffice,   Hence, 

the  iteration  !2,U,l)  with  the  fi. 's  chosen  according  to  the  criteria  of 

of  theorem  2„U,1  is  more  properly  classed  as  a  direct  method  than  as  an 

iterative  one,  and  direct  methods  are  net  the  principal  concern  of  this 

tnesis,   In  Chapter  3-  however,  we  will  use  tne  criteria  of  theorem  2, 4,1 

not  to  obtain  the  ft. 's  themselves,  but  to  obtain  certain  scalar 

acceleration  parameters  used  in  the  actual  iteration, 

2,5  Extrapolation  Matrices_When  the_Block  Matrices  Have  a  Common  Basis 
of  Eigenvectors 

Consider  now  the  linear  system  (2,1,1,  arising  from  the 
discretization  of  the  Dirichlet  problem  on  a  rectangular  domain  R  for  a 
self -aaj.omt ,  elliptic  partial  differential  equation  of  tne  special  form 
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^  f   x)  ^^zl)  _   d  (f  (y)  ^lls 

ox   1      Ox       Oy   2      oy 

-  (c, ,x:  -  a_(y))u(x5y)  =  g^x,y), 


(2.5.1) 


where  f  x   f  ;.  ,  a  (x),  o     y   are  continuous  in  the  closure  of  R 

and  satisfy  f  [x)  >  0,  f  (y)  >  0;  o  (x)  >  0,  o   (y)   >  0.   Let  a  uniform 

mesh  be  imposed  on  R,  and  suppose  that  the  mesh  points  are  numbered 

either  by  rows  or  by  columns  starting  from  the  upper  left.   Thus,  the 

mesh  points  .and  hence  the  unknowns  in  (2„l.l))  are  grouped  into  q  blocks 

of  r  elements  each.   It  follows  that  for  I  <  i  <  q,  the  blocks  B.  ,  n    and 

-   —  i,i-l 

B.  .  _  of  the  block  Jacobi  matrix  derived  from  A   are  square  of  order  r  and 

i  i+l 

all  have  a  common  basis  of  eigenvectors  u„e„,  a  set  of  r  linearly 

independent  eigenvectors.,,  In  this  case  the  matrices  [ft  }  can  be  chosen 

in  various  ways  so  that  they  share  this  common  basis  of  eigenvectors,  and 

consequently  a  useful  factorization  of  the  characteristic  polynomial  of 

£  is  possible o   In  Chapter  3  we  will  investigate  some  iterations  which 

are  equivalent  to  the  use  of  extrapolation  matrices  having  the  same 

eigenvectors  as  the  matrices  B.  .  .  and  B .  .  . .   The  matrices  {ft.}  will 

i.i-l      i,i+l  i 

then  be  determined  by  specifying  certain  of  their  eigenvalues. 

Let  t   .  .,  ,  6  .  .  , ,  and  ^  denote  the  eigenvalues  of  B.  .  _, 

i  .  i  - 1   i.i-l       i  i,i-l' 

B.  .  ,  and  ft.  respectively  to  which  the  I  r-dimensional)  eigenvector  i 
1 5 1+ 1      l 

belongs,  and  consider  the  scalar  iteration 

1         1,1-11-1      l,lrl  i     i' 

1  <  i  <  q„ 
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Let  Ju\    note  the  iteration  matrix  for  (2.5.2),  and  let  $  (\)  denote 

the  characteristic  polynomial  of  £r  ,  We  then  have  the  following 

decomposition  theorem,  to  be  used  in  Chapter  3° 

Theorem  2„5,,1;   If  the  matrices  B.  .  .. ,  B.  ,  _.  and  ft.  have  a  common 

<L L_^      1,1+1  1 

basis  of  eigenvectors  for  1  <  i  <  q,  then  the  characteristic  polynomial 

r 
of  £n  is  <|>(\)  =   n  <BJ(X). 

Proof  •   If  x  .:  [x  ,  x  ,  — ,  x  ]   is  an  eigenvector  of  iV\  belonging  to 
the  eigenvalue  \,  then  it  follows  from  (2,5 » 2)  that 

\x.  =  w?(\£.    x.  _  +  p.    -X.  _)  4   L-w?)x.,  1  <  i  <  q, 

1      11, 1-1  1-1      L , 1+1  1+1      %     11     -    - 

Hence  GJ  X  --■■   0,  where 


GJ  = 


[(1-^  -  . 


2  21 


q   q-q-1 


q-1  q-l,q 


L-ooJ)-\l 

q 


follows  that  A.  is  an  eigenvalue  of  JlTj  if  and  only  if  \  is  a  root  of 
0,  and  therefore  $J(\)  =  del 
Let  H  be  the  square  matrix  of  order  r  whose  j-th  column  is  %    , 

q 

and  let  H=  E©H,  where  ©  denotes  the  direct  sum.   The  linear 
independence  of  E  ,    £  ,  — -  £   implies  that  H  is  non-singular,  and  so  is 
letl  j)}   wn-re  G  was  defined  in  the  proof  of  theorem  2.^.1,  and 
it  follows  that  4>(\)  =  det  (H         G  is  block  tridiagonal ,  and  each 
non-null  block  of  G  is  square  of  order  r   Hence,  the  effect  of 


3h 


a  of  G  by  H  is  to  transform  each  non-null  block  of  G  by  H, 
and  diagonal! ze  it.   Thus, 


H~Vb, 


1,  l-l 


H  = 


and 


H  fi.B.    nH 


i   i ;  1-  1 


[(l-w})-\] 


H*  [(i-fi.)-X" 

[(l-<)-\] 

Hence,  H       is  an  shown  in  Figure  k,     . 

o   be  the  permutation  which  selects  the  first  element  of  the 
first  block  first,  the  first  element  of  the  second  block  second,  and  so  on 
through  the  first  elements  of  the  various  blocks^  then  likewise  through 
the  second  elements  of  the  various  blocks,  etc.   If  P  is  the  permutation 

_T   .Tl 

matrix  corresponding  to  a,  then  cj>'  ,}    ■■--   det(P  H   GHP   )„   But 
E  H         =  £  ©  GJ,  and  hence  d>(\)  =  IT  det(GJi  =  TT  0J(\)o 


j=l 


j-i 


Q,E0D, 
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3  c      SEQUENTIAL  EXTRAPOLATED  IMPLICIT  METHODS 


3,1     A  Model  Problem 

The   ususl  discretization  of  the  Dirichlet  problem  for  Poisson's 
equation 


dx  dy 


(3.1.1) 


on  a  rectangle  R  using  a  five  point  star  leads  to  a  linear  system  whose 
coefficient  matrix  A  is  block  tri-diagonal  with  tri-diagonal  diagonal 
blocks.   Let  a  uniform  r  x  q  mesh  be  imposed  on  the  interior  of  R,  and 
suppose  that  the  mesh  points  are  numbered  by  columns  starting  from  the 
upper  left „  Then 

D   -I 


■I   D 


where  I  is  the  identity  matrix  of  order  r,  D  is  of  order  v,    and 


h       -1 


D  ^ 


1   k 


The  resulting  linear  system  (2.1.1)  will  be  referred  to  in  what  follows 
as  the  model  problem,, 
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For  the  model  problem,  each  non-null  block  of  the  block  Jacobi 
matrix  derived  from  A  is  simply^  D  ,    and  the  block  Gauss-Seidel  iteration 
for  the  solution  of  (2.10l)  is 

m+1   ^-IrTT  ro+1  ^   m 


Z 


^x  =  D"X[Z  u+"  ,  Z  "   +  K.],   1<  i  <  q,  (3.1.2) 

i  i-I     i+1    i  '    -   -  *  w    ' 


or 


_J7  m+ 1   _  m+ 1   „  m     v  ,  /  _  ..  _  \ 

DZ.    =sZ._+Z._+K..   1  <  l  <  q.  (3.1.3) 

i       i-1     l+l    i'    -   -  * 

The  block  SOR  iteration  obtained  by  extrapolating  (3.1.2)  is  frequently 
called  SLOR  (successive  line  over-relaxation)  since  it  consists  of  updating 
simultaneously  the  unknowns  corresponding  to  an  entire  vertical  line  of 
mesh  poin~. 

3.2  The  Sequential  Extrapolated  Implicit  Method  (SEl) 

~*  m+1   ~*  m  ~~ 
If  Z     -  Z   =  Zj  the  exact  solution  of  the  model  problem, 

then  for  any  scalar  s, 

/      •  l~  m+1   ^"  m+ 1   ^m      Zm   rT        .  /     \ 

(D  +  sI;Z     =  Z  .  .  +  Z  .  _  +  sIZ.  +  K.  .   1  <  l  <  q,     (3.2.1) 


or 


Z  m+1  =  CD  +  si)"1  [Z  ^  +  Z  m  ,  +  sIZ  m  +  K. ].       (3.2.2) 
i  i-I     i+I      i    i  '  ' 


1  <  i  <  q. 
The  iteration  (3.2.2)  can  be  written  as 


z  m+1  =  qd"1  [z  mtl  +  z  m  n  h  it.  ]  +  (i  -  n)z  m        2.3) 

i  i-l     i+l    i  i  J       ~" 


i  <  i  <  0.  ; 
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where  ft  =  (D  +  si)"  D.   That  is,  (3=202)  is  the  VSOR  iteration  (2.1+.1) 


i  ■  "I 

with  ft.  =  (.D  +  si)   D„  1  <  i  <  q„   We  denote  the  iteration  matrix  for 

f3»2.2)  by  S. 

UJ 

The  most  efficient  way  to  perform  SLOR  is  to  use  (3.1.2)  to  compute 


1 
—  m  +  p 

Z  .    "  and  then  to  compute 


1 

2 
Z  1   =  <4Z.       "Zj+Z.^l^i^q,  (3=2,4) 


— .  ni  -t-  —   _ 
m+1   ,,rT     2   "mi   "m 


Hence,  the  additional  work  required  to  perform  SLOR  over  that  required 
for  Gauss-Seidel  iteration  is  one  multiplication  and  two  additions  per 
unknown, 

•ation  (3°2„2)  requires  that  (D  +  si)  "  be  computed  instead 
of  D   ,  but  both  of  these  matrix  inversions  require  the  same  amount  of 
work  since  the  diagonal  of  D  is  non-zero 0  Therefore  the  only  additional 
work  required  to  perform  SET  over  that  for  Gauss-Seidel  iteration  is  one 
multiplication  and  one  addition  per  unknown  in  order  to  compute  the  vector 
on  the  right  hand  side  of  (. 3o2„l)„   Hence  SEI  requires  one  less  addition 
per  unknown  per  iteration  than  SLOR„ 

Any  eigenvector  of  D  is  also  an  eigenvector  of  ft,,  and  since  the 
eigenvectors  of  D  form  a  basis,  the  converse  also  holds.   Hence,  the 
decomposition  theorem  205»1  applies. 

Let  d  be  an  eigenvalue  of  D  and  %    an  eigenvector  belonging  to 
d,   If  (3  denotes  the  corresponding  eigenvalue  of  D  ,  then  (3  =  —,  and  the 
eigenvalue  of  ft  =  (D  +•  si)   D  to  which  %   belongs  is  - — r~  ,  which  we 
denote  by  w(g;S)0 
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Let  JlT  denote  the  matrix  for  the  (scalar)  iteration 


x^+1  =  <o(p)p(x^  +  x*  )  +  (l-u(p))x»  L  <  i  <  q,     (3«2o5) 


and  let  ty  (\)  be  the  characteristic  polynomial  of  £T.  It  follows  from 
theorem  2.5.1  that  if  w(p)  =  w(p,s)  for  all  3,  then  the  characteristic 
polynomial  of  S^  is  ty(\)  =  TT\|r  (XJ,  where  the  product  is  taken  over  all 

eigenvalues  of  D 

ft 
Let  P-Vo'i  denote  the  maximum  of  the  moduli  of  the  zeroes  of 

ft  ft 

\|/  (\)  and  let  w  (ft)  denote  the  value  of  w(f3)  which  minimizes  p    fa\.      It 

follows  from  the  theory  of  SOR  that  for  w(p)  >  wh(P);  the  zeroes  of 

\|r  (\)  are  complex  of  modulus  (w(ft)-i),, 

For  fixed  s,  let  p(S,  ,)  denote  the   spectral  radius  of  S,  lC   If 

'CO  c  CO 

§_  and  ft  are  the  smallest  and  largest  eigenvalues  respectively  of  D  } 

then  p(S  )  =   max   __  p  ,        »,  and  we  let  p    =  min  p(S  ). 
R_  <  ft  <  ft     p'by  ox,x  s 

The  eigenvalues  of  D   are 


and  so 


P,  - 


U  -  2cos( 


rJIL^ 


r+l' 


\  <  &  <  Pj  <  P  <  | 


1  <  J  <  r  > 


Let 


A  = 


1   0 


be  square  of  order  q0   The  eigenvalues  of  A  are  CC.  =  2cos  —r  }    1  <  i  <  q_f 
1   since  the  block  Jacob!  matrix  derived  from  A  is  M  =  A  @  D   where 
@  denotes  the  tensor  product^  it  follows  that  the  eigenvalues  of  M  are 
\       a  a.B.j  1  <  i  <  q_s    1  <  j  <  r0   Moreover,,  the  theory  of  SOR  applied  to 
the  iteration  (3o2„5)  shows  that 


V-p)  ■  — 7—  .  (3.2.6) 


1  +  \[i-aap2 


where  OL   =a  max  {CH.}  „   Since  the  largest  positive  eigenvalue  of  M  is  ji  =  a  Pj 


we  have 


b    bVH; 


+\pV 


r:  VSOR  iteration  (2.U.1)  becomes  SOR  lif  fl.  a  w  I  1  <  i  <  q, 

Ls  case,,  theorem  205ol  shows  that  the  characteristic  polynomial  of 

8  •  B 

£0      is  <£>(a)  3  TI  O^'.a.1  j,  where  cj>  (A.)  is  the  characteristic  polynomial 

a      b  B  o 

of  JlT  with      1  to    The  maximum  of  the  moduli  of  the  zeroes  of  $  (\) 


B 

J 

t 


is  pw  , 


Let   s^tO    =:  -  b~7pT  -   x   ?    and  denote   sb(p)by  sbD      Then  for 


Lemma 

j_b 
B  <  B-  ub(p)  <  w(B,sb)  <  wfe(3)  a  u^, 

oofs   s  (6)  <  0  since  ^.(B)  >  1,  and  so  1  +  s  J  <  1  +  s  B  if  0  <  B  <  P, 

Hence,  w(B;s  )  <  wb(B)0   Substituting  (3„2„6)  into  the  formula  for  s  (p) 

and  differentiating  shows  that  -rr   s.  (b)  <  0.  and  hence  s,  <  s,  (b)  for 

dB  b       '  b    b 

B  <  B .   Then 


b  hs 


an 


d  so  u)  (B)  <  w(B5s,  ) 


Q.E.D, 


Ul 


Theorem   3.2.1:      If  s  =   s,  ,    then: 
' b 


R  R  R 

(i)  r  (X.)  =  $  (A.),  and  hence  p^  ^  =  pg0R 

b ' 

(it)     ForP<e,   p^(p)  <ol{&)%)  <PS0R 

(iiij   lim  p  ,    \  =  0  for  any  s. 
R.*0   ^>S) 


Proof : 


(i)  By  definition,  s  is  that  value  of  s  such  that 

—        —  —  6*      R 

<*>(0,s)  =  wb(P)»   But  wb(P)  =  °*h>    and  so  V  (M  =  $  (^).   Since  the  moduli 

o  

of  the  roots  of  <$>   (\)    are  equal  to  pC(~D  for  all  B  <  (3,  it  follows  that 

B 

Pw  (0)  -  PSOR: 
b 

(ii)   It  follows  from  lemma  3«2.1  and  the  theory  of  SOR  that  the 

roots  of  ty  (\)  are  complex  of  modulus  (a)(f3,s  )-l)  for  0  <  R,  and 

(w  (B)-l))  <  (w(0,s  )-l)  <  (w^P)-l).   The  assertion  (ii)  then  follows 

from  the  definitions  of  the  quantities  involved, 

(iii)   lim  w(B,s)  =  1,  and  since  p  ,    \  =  (w(0,s)-l),  we  have 
fi->0  \P>S. 

lim  p  /  „\  =  0, 

B^O  W(^S)  Q.E.D. 


Corrollary  3,2,1:    PSEI  *  PS(}R  . 

Proof;   The  theory  of  SOR  applied  to  J^  shows  that  P^/q   \  is  minimal  if 

and  only  if  w(r  s)  =  w  (r)  which  is  the  case  if  and  only  if  s  a  s.  ,   It 

b  b 

then  follows  from  part  (ii)  of  theorem  3.2.1  that    max    Pm(a   1  ^s 

R_  <  R  <  R    lP'Sj 
minimal  if  and  only  if  s  =  s  ,  and  part  (i)  of  the  theorem  shows  that 


PSEI  ''"'   pS0R 


Q.E,D, 


1*2 

t  |  '  normalized  to  have  (Euclidean)  length  1  be  the 

-1  ~"*  I  r 

eigenvector  of  D   belonging  to  P.j  1  <  J  £  r°   ^he  vectors  { |   }   form 

J  -1- 

an  ortho-normal  basis  for  r-space  since  D   is  symmetric.,  and  indeed 

'r+1  '     r+1  '    3  r+1 

t  e  "  denote  the  q-dimensional  column  vector  having  1  in  the  i-th 
position  and  01  s  elsewhere,,   The  vectors  [e   }  _r  form  an  ortho-normal  basis 
for  q-space,,  and  therefore  the  vectors  [e  '  @  £  '\)^-S.^-S.(:lj^-<S.^'S.x's 
form  an  ortho-normal  basis  for  qr-space,,  where  ©  denotes  the  tensor 
product  o 

eigenvector  of  D  "  corresponding  to  f3  is  |   ,  and  it 
follows  from  part  (i)  of  theorem  3.2,1  that  if  the  initial  error  E   lies 
in  the  subspace  spanned  by  {e    @  %      }      then  the  error  E   in  the  m-th 
iterate  for  SE.r  is  identical  to  the  error  in  the  m-th  iterate  for  SLOR,,  Vm, 

Because  of  part      of  the  same  theorem,  however,,  it  is  to  be  expected 

-*  0 
that  if  E   does  not  lie  entirely  in  this  subspace,,  then  a  fixed  number  of 

iterations  with  SEI  will  produce  a  greater  reduction  in  the  I     norm  of  the 

error  than  the  same  number  of  iterations  with  SLOR,,   Numerical  results 

which  support  this  expectation  are  presented  in  Section  3-^. 


3 . 3  Cv clic  Chebyshev  SEI 

The  cyclic  Chebyshev  semi -iterative  method   [6  j  11,  Chap „  5  ]  for 
the  solution  of  the  model  problem  is 

z.   =  u)    d  [z._+z,_+K.J  +  (l-w   ;z  . ,  (3,3,1a; 
1  1-1     1+1    1  1 

1  - :  t,  3 ,,  2 ,,     '  J 


><3 

Zm+1  =  c2m+2  D4[Zm+!  +  Z  "*}  +  K.]   +    (l-^m+2)Zm,         (3.3.1b) 
1  l-l  l+l  i  l' 

i  =  2,    h,    6,    ---  , 


r     mi 

where  the  sequence  [u   J  satisfies 


1   ,    2     1 

W   =  1.    W   = 


) 


-2 

1-  L 


1-g-eJ" 


(3.3.2) 


and  u  =  p(M)  , 

In  what  follows  (3»3.l)  with  the  sequence  {w  }  satisfying 


(3.3.2)  will  be  called  cyclic  Chebyshev  SLOR,  while  (3.3.l)  with 


m 

10   =  OJ 


b' 


Vm,  will  be  called  SLOR  with  the  odd-even  ordering.   Iteration  (2,2.3) 

with  w  =  to  will  be  called  SLOR  with  the  natural  ordering,, 
b  — " 

If  the  sequence  fw  }  satisfies  (3.3.2),  then  lim  w  =  w   and 
^      l  '  \ -><■->      I?  m->°°      b 

it  follows  that  the  asymptotic  rate  of  convergence  of  cyclic  Chebyshev  SLOR 
is  equal  to  that  of  SL0R<   However,  if  E   denotes  the  error  in  the  m-th 
iterate  for  cyclic  Chebyshev  SLOR  and  if  E   denotes  the  error  in  the  m-th 
iterate  for  SLOR  (with  any  consistent  ordering),  then  Golub  and  Varga  [6] 
have  shown  that 


max  | |E  m| |  <  max| |E  m| | ,   m  <  2, 

;o        -  o 
z         z 


where         denotes  the  I     norm.  Moreover,  the  sequence  (||E   ||)  is 


mi 


strictly  decreasing.,  whereas  the  sequence  (j|E    j}  may  increase 
initially. 


tsm}  by 


Let  the  sequence  {w  )  satisfy  (3=3=2),  and  define  a  sequence 


Sm  =  I   [-i  -  1],   m>  1  .  (3,3,3) 


For  0  <  &  <  f>.  define  the  sequence  {w  (f3,s  )}  by 

m,„   mN      1 


>  (p,sm)  -    xm  ,  m>  l.         (3.3.U) 

1  +  si 


Then 


m/o  „m\    m  v_/ 


The  iteration 


m- 1    /  2m-lM-l  r   m  _  m      2m+l  _  m  v    -, 

Z       I D  +  s    I J    Z  ,  n  +  Z  ,  ,  +  s     Z  .  +  K  ,  J , 

i  l-l  l+l          i  i  ' 

1  =  13  5—  (3<3'5a) 

~Z  m- 1        /_  2m+27>-l   rr"  m+1  ^  nul         2m+2  r"  m  IT  1 

Z                   DtS          1)        [Z-t+Z        -,+s  Z    .    +  K,  J 

i             v  l-l  l+l                          l  i 

i    .  2,   k,  6,   -  (3.3.5b) 


for  the  solution  of  the  model  problem  will  be  called  cyclic  Chebyshev  SEI, 
Iteration  ( 3  -  3  -  5 )  with  s   replaced  by  s,,  Vm.,  will  be  called  SEI  with  the 
odd -even  ordering;  while  (3,2,2)  with  s  -  s  will  be  called  SEI  with  the 
natural  ordering. 

The  remarks  in  Section  3  2  concerning  the  convergence  properties 
of  SEI  relative  to  those  of  SLOR  apply  also  to  the  convergence  properties 
of  cyclic  Chebysnev  SEI  relative  to  those  of  cyclic  Chebyshev  SLORD 


^5 

3.^4       Leal  Results 

Numerical  experiments  were  performed  to  compare  the  actual 
rates  of  convergence  of  SLOR  and  SEI  with  both  the  natural  and  I 
odd-even  ordering,  cyclic  Chebyshev  SLOR,  and  cyclic  Chebyshev  SEI.   The 
test  problem  selected  was  the  model  problem  with  f(x,y)  =  0  (i.e., 
Laplace's  equation)  and  homogeneous  boundary  conditions.   The  rectangle 
R  was  taken  to  be  a  square  containing  N  mesh  points  on  a  side.   For 
N  =  10,  20,  30,  and  Uo  the  number  of  iterations  with  each  method  needed 
to  reduce  the  I     norm  of  the  error  below  1  was  determined. 


Starting  vectors  were  constructed  as  follows.   Let  Q,  be  the 

N-dimensional  vector  whose  j-th  component  is  (-l)     ,  and  let  J  be 

_    N   _  , 
the  N-dimensional  vector  whose  j-th  component  is  j.   Let  r\  -     ^    £  , 

^1^2       -*  N  -1        ^=1 

where   |  ,  £  s    ,  £   are  the  eigenvectors  of  D   .   Experiments  were 

performed  using  each  of  the  three  starting  vectors  Z   -  J  ©  Q  , 

Z   =  10  Q  ©  Q„  and  Z   =  J  ©  r\. 

The  graphs  shown  in  Figures  5 ,    6,  and  7  were  constructed  by 

interpolating  linearly  between  the  data  points,   It  can  be  seen  from 

these  graphs  that  SEI  converges  much  more  rapidly  than  SLOR  for  certain 

starting  vectors  but  that  the  relative  advantage  of  SEI  over  SLOR  is 

strongly  dependent  on  the  starting  vector  chosen. 
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3 . 5  The  Variable  Sequential  Extrapolated  Implicit  Method  (VSEl) 

Consider  the  following  generalization  of  (3.202)  for  the 
solution  of  the  model  problem. 


_  m+1   , ,  ,_     T\-l  r_  m+1   _  m     T.      „  ran 
Z.     u)  (D  +  s.I  )    [Z,_+ZJ_+K  +  s.Z  .  J 

1      i      l        l-l     i+l    i    li 


+  (l  -  u.)  Z  .  ,   1  <  i  <  q; 


(3.5.1) 


where  (w.),  and  {s.},  are  any  two  sequences  of  real  numbers,   The  special 
case  of  (3.5.l)  in  which  w_  =  ]_,,  Vi,  will  be  called  1-parameter  VSEI 
whereas  the  general  case  will  be  called  2-parameter  VSEI „ 
Iteration  (3.5=l)  can  be  written  as 

I*  m+1     n-lrr*m+l   I'm  -,    /_     \^  m 
Z  .     u.u        IZ._+Z._J+(I-fi.;Z.J 

1      1       l-l     i+l  1   1 

(3.5.2) 
1  <  i  <  <L, 

where  SI     =   w.(D  +  s.l)   D„   Hence  ( 3.5.1*)  is  a  special  case  of  the  VSOR 
111 

iteration  (2,U.l)0   The  iteration  matrix  for  (3«5°l)  will  be  denoted  by 

The  work  required  to  perform  1-parameter  VSEI  is  the  same  as 
for  SET,   The  most  efficient  way  to  perform  2-parameter  VSEI  is  to 
compute 

1 

p»        2        fT.  _  \  - 1    r  _  m+ 1        _   m  _   m       „-, 

Z  =   '  D  +   s.  I)        [Z        n+Z.^+sZ.4-K.J, 

1  i-1  i+l  11  1 


(3.5.3) 


1  <  i   <  q   , 
and  then  to   compute 


SO 


1 

z  '"  =  ".  [Z  .m+2  -  Zm]  +  z  m,  l  <:  i  <  q  ,         (3.5.U) 

ii      i     i'   —   — 

e   the  additional  work  required  to  perform  2-parameter  VSEI  over  that 

required  for  SLOR  is  one  multiplication  and  one  addition  per  unknown  per 

iteration. 

The  eigenvectors  of  W.(D  +  s.l)   D  are  the  same  as  those  of 
°  ii 

-1  - 

D  j,  and  so  the  decomposition  theorem  2„5ol  applies.   If  £  is  an  eigen- 
vector of  D  "  belonging  to  the  eigenvalue  B,  then  the  eigenvalue  of 
D  +  s.l)   D  to  which  f  belongs  is 

i 


1  -*  s.B 


which  we  denote  by  w . (B,s.) 


Let  £L  denote  the  matrix  for  the  (scalar)  iteration 
ft 

xm+1  =  uWB,s.)e(Xm^  f  xm  J  4  (l-U).(p,s.))Xm, 

1      11    1-1    1+1'       11   1 


i  < .  i  <  q.  , 


(3.5.5) 


R  R 

and  let  a-  (\)  denote  the  characteristic  polynomial  of  £  „  It  follows 

from  theorem  2,5=1  that  the  characteristic  polynomial  of  S   is 

R  -1 

•  IT  ^  (A-)?  where  the  product  is  taken  over  all  eigenvalues  of  D 


Let  p0  denote  the  spectral  radius  of  Jl~  „  and  let  PvqT?-  denote 

the  spectral  radius  of  Sn„  Then  p. .__...  =   max      p  ,  where  8  and  B 

£  <  3  <  F 
denote  the  smallest  and  largest  eigenvalues  of  D  "  respectively , 

w  (B,s  )  is  specified  for  two  distinct  values  of  8,  8n  and 
l    l  1 

8„  say,  then  s„  and  w  are  determined  by  the  formulas 
2  ii 
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w.  :  (l  4  siP1)  uii'S1,s.)  .  (3,5.6b) 


P     6 
3,6  VSEI  with  d^1  -  pQ2    --   0 


The  criteria  of  theorem  2cUci  with  fi  =  w/ft  s  )  and 

1      L   '  1 


q 


B.  .  n  =  B.    ,  =  P  can  be  used  to  determine  sequences  fw.(Pn3s  )}_,  and 
i,i-l    i,i+l  ix  1J  i  1 

-  q  Pl      P2 

&  ,s.,').,  such  that  £   and  £,   are  nilpotent;  icec1  such  that 

^        J-        -L  06  Ob 


to  be  used  in  [3.5. l)«   It  then  follows  from  (.3.5  »6b)  and  theorem  2C4,1 


o   =  p   ---  0o   Then  (3»5=6)  determines  the  sequences  {s.}   and  l.^."1 


that  I -parameter  VSEI  is  the  special  case  (3.  =  0  of  2-parameter  VSEI0 

It  can  be  shown  that  p   is  the  same  function  of  p  if  the 

criteria  of  either  part  (i)  or  part   11,,  of  theorem  2,4,1  are  used  to 

determine  (w.(P,ss.)l  and  fw.(P  ,s.)}  .  in  which  case  (3c5„l)  will  be 
i  '  1"  £  '1       i   2  i   1 

called  uni -directional  VSEI.   If  the  criteria  of  part  (iii)  of  the 
tneorem  are  used,  (3.5.l)  will  be  called  bi-directional  VSE[ 

/.'as  found  experimentally  that  if  P   :  P   then   max    p 

increases  as  P   increases  while    max      p   decreases,   Moreover 

PL  <  P  <  P    Q 

P    P 
max     p„  -:  pE-.   it  follows  that  if  p..  is  the  values  of  p.  such  that 

P  <  6      "  "  lb  X 

_    _         _      ^ 

P  P  P 

max  p     =  pE;,    then       max  p^   is   minimal  for  Pn    -  pni  , 

Px  <  P  <  P     °  °  &  <  P  <  P  X  lb 


For  the  model  problem,  p    =  p^  /-v,  where  p^/R\  was  defined 

b 

3      3 


nor 


in  Section  3=20   The  theory  of  SOR  shows  that  p   ,->      ,— \ 

-dp 

0  <  P  <  Pj  and  —  p^  ,r\  =  +  oo  for  p  =  P  +  „   The  approximate  shape  of 

the  graph  of  p  as  a  function  of  P  is  shown  in  Figure  8,  and  the  graph  of 

bh  (3  =  Po  is  shown  for  compari! 
b 


Fj  gure  8 


For  the  model  problem.  §  >  ?-  and  (3  <  —j  Vr„   Moreover 


lim  £  =  g-  and   lim  p   -    Hence  PYSEJ 


max 


r 


r  -■>«. 


ft.  <  P   <  & 


u 


P  P 

p     <  max          p    . 

~  1  P        -1 

6  -  p 


iteration    I  s)    with  u.(3,s. )    =  ; ■ — ,    1  <  l  <   q.    was   performed    to 

i  '  l    1  +  s.P;    ~   - 

! 

Q 

determine  p   as  a  function  of  P  for  q  =  10,  19,  and  33,   Formula  (3.5.6) 

with  P0  =  ~   W5£-  used  to  determine  ( s . }  _r  and  fw.}  ,   For  the  ^-parameter 
c.  1  J       i  J 

case  p   was  determined  by  numerical  search.   The  results  of  these 


!■■■  r  '  m<  tit  J:  •  x  e  ih  esented  Ln  Dab]  <  1-U 
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Table  1.   1-Parameter,  Uni-Direction  VSEI 


» 

PVSEI 

PS0R 

^SEI 

RS0R 

10 

.55 

.56 

.60 

.58 

19 

.7»» 

.73 

,30 

,31 

33 

-85 

=  83 

.16 

.19 

Table  2,   1-Parameter,  Bi-Directional  VSEI 


q 

PVSEI 

pS0R 

RSEI 

RS0R 

10 

M 

,56 

.82 

=  58 

19 

,65 

=  73 

M 

=  31 

33 

o78 

,83 

=  25 

=  19 

Table  3,   2-Psrameter,  Uni-Directional  VSEI 


q 

\b 

PVSEI 

PS0R 

RVSEI 

RS0R 

10 

.Ul 

M 

.,56 

.81< 

-58 

19 

A7 

,6k 

.73 

M 

.31 

33 

.*9 

=  78 

=  83 

=  25 

,19 
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Table  hc      2«Parametern  Bi-Directional  VSEI 


q 

pib 



PVSEI 

pS0R 

^SEI 

r        i 

SOR   i 

10 

.3^ 

o32 

.56 

1.1 

,a 

19 

.hi 

s   .  .... , 

,73 

o62 

•31 

1  r> 

M 

•     ,69     .83 

1 .i .... 

.37 

=  19 

For  bi-directional  VSEI,,  best  results  were  obtained  with  m  in 
part   1  i   of  theorem  204ol  equal  to  the  greatest  integer  not  exceeding 
'^r~i  and  the  results  presented  are  for  this  case.   In  all  cases  except 
that  of  1-parameter,  uni-directional  VSEI,  the  rate  of  convergence  of 
VSE]  for  the  model  problem  exceeds  that  of  SLOR,  and  in  the  case  of 
2-parameter  VSE;r  ,  the  improvement  is  substantial. 


.   CONCLUSION 

Summary 

Several  schemes  have  been  investigated  for  impn      the  rate 
of  convergence  of  extrapolated  Gauss-Seidel  iteration  by  the  introduction 
of  a  multiplicity  of  extrapolation  parameters  to  replace  the  single 
scalar  parameter  used  by  SOR.   In  Chapter  2  it  was  shown  that  the  use  of 
two  optimally  chosen  extrapolation  factors  results  in  an  improved  rate 
of  convergence  for  linear  systems  whose  coefficient  matrices  are 
consistently  ordered  S-matrices0   In  the  case  of  linear  systems  arising 
in  the  numerical  solution  of  boundary  value  problems  for  elliptic- 
partial  differential  equations  the  improvement  is  small  because  u_  for 
these  systems  is  small.   It  is  evident,  however,  that  there  exist  linear 
systems  for  which  u_  and  u  are  more  nearly  equal,  and  for  such  systems 
the  use  of  two  optimally  chosen  factors  offers  a  substantial  advantage. 

It  was  also  shown  in  Chapter  2  that  by  use  of  matrices 
rather  than  scalars  as  extrapolation  parameters,  an  extrapolated 
Gauss-Seidel  iteration  having  an  infinite  rate  of  convergence  can  be 
constructed,  although  its  implementation  requires  more  work  per 
iteration  than  extrapolation  by  a  scalar. 

The  SEI  and  VSEI  methods  of  Chapter  3  were  analyzed  for  a 
more  limited  class  of  linear  systems,  namely  those  arising  from  the 
discretization  of  the  Dirichlet   problem  for  (3.1.].)  on  a  rectangu] 
domain ,      In  this  case  SEI  was  shown  to  be  equivalent  to  several 
simultaneous  SOR  iterations,  each  on  a  different  sub 

a  different  scalar  extrapolation  factor.  Moreover,  SEI  was  shown  to 
have  convergence  properties  superior  to  those  of  SOR  for  the  problems 
considered  and  to  require  less  work. 
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.-  investigation  of  VSEI  was  largely  experimental,  but  the 
results  indicate  that  for  the  class  of  problems  considered  it  is  superior 
to  SOR  not  only  with  respect  to  average  rate  of  convergence  but  also  with 
respect  to  asymptotic  rate  of  convergence. 

kc2     Extensions 

It  seems  probable  that  the  results  of  Chapter  3  can  be 
extended  to  linear  systems  associated  with  the  more  general  self-adjoint,, 
elliptic  partial  differential  equation  (2„5.l)  since  the  non-null  blocks 
of  the  Jacobi  matrix  in  this  case  still  have  a  common  basis  of 
eigenvectors  though  different  eigenvalues.   Extensions  to  boundary  value 
problems  for  non-rectangular  regions  appear  more  difficult,  but  such 
extensions  would  be  very  desirable. 
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